#clean up
rm(list=ls())

#options
options(scipen=8)

#Number of samples
N<-100

#### PROCESSING FORECASTS ####
for(i in 1:N){
	check<-file.exists(paste('bootstrap/forecast',i,'.csv',sep=''))
	if(check){
		assign(paste('f',i,sep=''),read.csv(paste('bootstrap/forecast',i,'.csv',sep='')))
	}else{
		print(paste('File forecast',i,'.csv does not exist.',sep=''))
		}
	}

#merge and drop component files
forecasts<-rbind(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,f25,f26,f27,f28,f29,f30,f31,f32,f33,f34,f35,f36,f37,f38,f39,f40,f41,f42,f43,f44,f45,f46,f47,f48,f49,f50,f51,f52,f53,f54,f55,f56,f57,f58,f59,f60,f61,f62,f63,f64,f65,f66,f67,f68,f69,f70,f71,f72,f73,f74,f75,f76,f77,f78,f79,f80,f81,f82,f83,f84,f85,f86,f87,f88,f89,f90,f91,f92,f93,f94,f95,f96,f97,f98,f99,f100)
rm(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,f25,f26,f27,f28,f29,f30,f31,f32,f33,f34,f35,f36,f37,f38,f39,f40,f41,f42,f43,f44,f45,f46,f47,f48,f49,f50,f51,f52,f53,f54,f55,f56,f57,f58,f59,f60,f61,f62,f63,f64,f65,f66,f67,f68,f69,f70,f71,f72,f73,f74,f75,f76,f77,f78,f79,f80,f81,f82,f83,f84,f85,f86,f87,f88,f89,f90,f91,f92,f93,f94,f95,f96,f97,f98,f99,f100)
forecasts<-subset(forecasts,subset=predictive.mean>=0 & predictive.mean<=100)
forecasts<-subset(forecasts,subset=STATE%in%c("Alaska","Hawaii"))

#Obtain sample for making predictions
subset.data<-read.csv('bootstrap/train1.csv')

#coefficient order
lm(subset.data$ideology~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
		subset.data$inc14+subset.data$cathOrth+subset.data$consRelig+subset.data$musJew+
		subset.data$mainline+subset.data$rural+subset.data$ownership+as.factor(subset.data$empstat)+
		subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2))
		
#make a forecast
cor.ak<-rep(NA,40)
cor.hi<-rep(NA,40)
se.ak<-rep(NA,40)
se.hi<-rep(NA,40)
for(i in 1:40){
local.northings<-forecasts$northings-450+50*(i-1)
kriges<-krige.conv(coords=cbind(subset.data$eastings,subset.data$northings), data=subset.data$ideology, locations=cbind(forecasts$eastings,local.northings),
			krige=krige.control(type.krige='sk', trend.d=~subset.data$age*subset.data$educ+as.factor(subset.data$race)*subset.data$female+
		subset.data$inc14+subset.data$cathOrth+subset.data$consRelig+subset.data$musJew+
		subset.data$mainline+subset.data$rural+subset.data$ownership+as.factor(subset.data$empstat)+
		subset.data$eastings*subset.data$northings+I(subset.data$eastings^2)+I(subset.data$northings^2),
			trend.l=~forecasts$age*forecasts$educ+as.factor(forecasts$race)*forecasts$female+
		forecasts$inc14+forecasts$cathOrth+forecasts$consRelig+forecasts$musJew+
		forecasts$main+forecasts$rural+forecasts$ownership+as.factor(forecasts$empstat)+
		forecasts$eastings*forecasts$northings+I(forecasts$eastings^2)+I(forecasts$northings^2),
			beta=c(45.6968, 0.1866,-1.4007,-9.0599, 0.5211,-4.0114, 0.2099, 6.9414, 16.3552,-5.0293, 4.8479, 0.8307, 4.9430,-1.3826,-0.2340, 0.2705/1000,-2.1224/1000,-0.6419/(1000^2),-1.3786/(1000^2),-0.0269, 6.7673,-2.1495,-0.6522/(1000^2)), 
			cov.pars=c(87.0752,4860.5613), 
			nugget=I(610.5291), cov.model='gaussian'))
ak.mod<-lm(scale(forecasts$predictive.mean[forecasts$STATE=="Alaska"])~scale(kriges$predict[forecasts$STATE=="Alaska"])-1)
hi.mod<-lm(scale(forecasts$predictive.mean[forecasts$STATE=="Hawaii"])~scale(kriges$predict[forecasts$STATE=="Hawaii"])-1)
cor.ak[i]<-ak.mod$coefficients
se.ak[i]<-vcov(ak.mod)
cor.hi[i]<-hi.mod$coefficients
se.hi[i]<-vcov(hi.mod)
}
min(cor.ak);max(cor.ak);cor.ak[10]
se.ak[which.min(cor.ak)];se.ak[which.max(cor.ak)];se.ak[10]
min(cor.hi);max(cor.hi);cor.hi[10]
se.hi[which.min(cor.hi)];se.hi[which.max(cor.hi)];se.hi[10]
